#---------------------------------------------------------------------------------------------
# Replication script for replicating results on competition between lobby clients
# In Egerod and Junk (2021)
# 
# The output of this script will correspond to Figure 4 in the article
# 
#--------------------------------------------------------------------------------------------

# set working directory to script location
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# load packages
library(readr);library(dplyr); library(ggplot2)
library(haven)
library(sp); library(spdep);  library(splm); library(igraph)

# load data
tot <- read_csv("LobbyFirmRevenue_May2020.csv")
w_mat <- read_csv("LobbyFirmWeightsMatrix.csv")
w_mat <- mat2listw(as.matrix(w_mat), style = "W")
tot$sd_btwn <- scale(tot$btwn) # standardize

# estimate model
slx.biv <- lagsarlm(log(tot.x+.5) ~ log(clients+.5)+ sd_btwn + 
                      log(tot.y+.5) ,
                    data = tot, w_mat,
                    method="eigen",  zero.policy=TRUE)
summary(slx.biv)

# bootstrap to get indirect/direct effects
set.seed(182)
lob_fx <- impacts(slx.biv, listw = w_mat, zero.policy = T,
                  R = 500, zstats = T, robust = T,
                  type = "HC1")
summary(lob_fx)

# some quantities of interest
slx.biv$rho
w.slx.biv$rho

# a regression table
stargazer(slx.biv,
          covariate.labels = c("Total # Ties (logged)",
                               "Betweenness Centrality",
                               "2014 Revenue (logged)"),
          title = "Linear SAR Coefficients:  Lobbying Firms Working for the Same Client",
          dep.var.labels = "2015 Revenue (logged)",
          add.lines = list(c("Spatial Parameter", round(slx.biv$rho, digits = 3)),
                           c(" ", round(slx.biv$rho.se, digits = 3))),
          out = "lobby_firm_regressions.html",
          type = "html")

#prepare for plotting

lob_ind <- t(apply(lob_fx$sres$indirect,2,
                   quantile,probs = c(0.025, 0.05, 0.5,
                                      0.95, 0.975)))

lob_dir <- t(apply(lob_fx$sres$direct,2,
                 quantile,probs = c(0.025, 0.05, 0.5,
                                    0.95, 0.975)))

lob_sims <- as.data.frame(rbind(lob_dir,lob_ind))
names(lob_sims) <- c("lwr95", "lwr90", "pe", "upr90", "upr95")


lob_sims$type <- c(rep("A: Direct Effect",3),
                   rep("B: Indirect Effect",3))

lob_sims$xvar <- c("connect", "btwm", "rev14") 
lob_sims$ordering <- 3:1


ann_text <- data.frame(pe = c(0.3,0), xvar = lob_sims$xvar[1:2], ordering = c(2.09,2.09), 
                       lab = c(round(slx.biv$rho,3),NA), lab.se = c(round(slx.biv$rho.se,3),NA),
                       type = c(lob_sims$type[1],lob_sims$type[4]))

ann_text2 <-data.frame(pe = c(0.3,0), xvar = lob_sims$xvar[1:2], ordering = c(2.05,2.05), 
                       lab = c(round(slx.biv$rho,3),NA), lab.se = c(round(slx.biv$rho.se,3),NA),
                       type = c(lob_sims$type[1],lob_sims$type[4]))

ann_text$cors <- paste(c("rho==", NA), c(ann_text$lab),
                       sep = "")
ann_text$cors[c(2)] <- NA
#ann_text$lab.se <- format(ann_text$lab.se, scientific=F)


ann_text2$SEs <- paste("(", ann_text$lab.se, ")",
                       sep = "")
ann_text2$SEs[c(2)] <- NA

# create the plot

p_client <- lob_sims%>%filter(xvar != "connect")%>%
ggplot(aes(x = pe, y = ordering)) +
  geom_point() +
  geom_errorbarh(aes(xmin = lwr95, xmax = upr95), height = 0.05,
                 colour = "black", alpha = 0.3) +
  geom_errorbarh(aes(xmin = lwr90, xmax = upr90), height = 0,
                 colour = "black") +
  theme_bw() + 
  geom_vline(xintercept = 0, lty = 3) +
  labs(x = "Estimate", y = NULL) +
  scale_y_continuous(limits = c(0.75, 2.25),
                     breaks = c(1,2),
                     labels = c("2014 Revenue","Betweenness")) +
  geom_text(data = ann_text, aes(label = cors), size = 3, parse = T) +
  geom_text(data = ann_text2, aes(label = SEs), size = 3) +
  facet_wrap( ~ type, scales = "free_x") 

p_client

ggsave(plot = p_client,
       filename = "LobbyClientResults.jpg", device = "jpg",
       width = 6.99, height = 4.99) 





